Theoretical and semi-analytical simulation for a two-predator-one-prey model during the mating period

The article introduces a new application which is a system of equations of two predators and one prey with the term of interaction between male and female of predators and prey. Such term appears when male and female of predators feed on the same prey during their mating period. The mathematical model has been studied theoretically and semi-analytically. The positivity, boundedness, local and global stability are proved for the system. The logarithm of multistage differential transform method (MsDTM) is used to study this new application. The MsDTM is used because it globally converges to the solution, it is a highly accurate, fast and simple approach. The stability analysis as well as semi-analytical solutions of the system are obtained to understand the dynamic of the model. Moreover, the effects of several parameters in the system are presented. As a results, we obtain the periodic solution when when the growth rate of prey is larger than the growth rate of both type of predators.


Introduction
Mathematical models phrased in terms of differential equations are very useful tools to describe several problems and applications in science such as biology [1], physics [2], fluid dynamic [3], medical science [4,5], etc. They are helpful to control the spread of infection, to understand the behavior of population, to observe the effect of several factors in population or to predict the solutions for different problems in the world. Some populations are described by a system of ordinary differential equations (ODEs) that introduces the interaction of different species [6].
The predator-prey model is a well known mathematical model and it was first introduced by Lotka [7] and Volterra [8]. It has been used widely to study the dynamic of interaction of two species in population ecology and associated ecological interactions [9]. In addition, researchers modified the model to study different biological and ecological scenarios. For example, they study the model among vertebrates and invertebrates in freshwater systems associated with fisheries [10] and with prey refuge and delay [11]. The modified models have been utilized for studying the impulsive effects [12], the impact of prey herd dread from the two types of predators [13], effect of food sharing among predators [14], the effect of two different food chains [15], studying the stationary probability distributions of its population densities [16], studying with disease in prey [17] and so on.

PLOS ONE
PLOS ONE | https://doi.org/10.1371/journal.pone.0289410 August 9, 2023 1 / 11 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 One of the behaviors of some type of animals such as the discus fish is that the male stays with the female during mating period to protect her. During this period both male and female feed together in same prey. The model that describes the interaction between one predator and two preys has been studied in [18,19], while the model of two predators and one prey has been studied in a few papers such as [20,21]. Al Qudah modified the predator-prey model by considering the effect of the term XY 2 , where Y 2 denotes the intersection between a male and a female of predators who feed on same prey (X) during mating period [22]. Then, Aljahdaly and Al qudah found the analytical solutions and discussed the biological relevance of the solutions [23]. The numerical solution was found by the exponential time differencing method [24]. However, no contributions in the existing literature have studied the effect of the term XY 2 on the two-predators-one-prey model.
Due to the difficulty of computing exact solutions for most ODEs that described real applications, semi-analytical solutions are found instead of exact solutions in explicit form [25]. Our previous work [26] presented the predator-prey with the effect of the term XY 2 and computed the semi-analytical solution by the differential transform method (DTM) [27], Adomian decomposition method (ADM) [28], 4 th Runge-Kutta (4 th RK) [29] and the multistage differential transform method (MsDTM). We proved that the MsDTM is more accurate than DTM and ADM and faster than the 4 th RK. In fact, semi-analytical methods find the solution in the form of series and the other terms of series are obtained from given initial conditions. Some semi-analytical techniques are based on Taylor expansion about the initial point t 0 such as Legendre Pseudospectral Scheme [30], residual power series method (RPSM) Laplace decomposition method (LDM) [31], DTM and ADM. These methods enable to find the solution in a very small time domain because the Taylor expansion is locally convergent. In other words, the DTM can find the solution only in small neighborhood about the initial point and can not obtain the solution in long domain. The DTM has been improved by using Padé approximation and by applying the multistage technique. Both improvements are able to find the solutions in long domain. However, the MsDTM is a reliable and optimal method to approximate solutions in a long domain for linear and nonlinear ODEs [26]. The MsDTM is a modification of the DTM by dividing the original time domain into sub-domains and applying the DTM method at each subdomain. In addition, the nonlinearity is treated by transforming the nonlinear term to a special summation based on differential transformation form. The accuracy of the MsDTM is controlled by the iteration number and time step size. The accuracy for the 4 th RK is controlled only by the time step size, therefore, the MsDTM can be faster than 4 th RK [32].
The novelty of this paper is presenting a two-predator-one-prey mathematical model during the predator mating period. In addition the paper studies the model theoretically including a positivity, boundedness, as well as local and global stability. The paper introduces the MsDTM logarithm and interprets the results.
The paper is organized as follows: section (2) proposes the mathematical model, section (3) shows the positivity and boundedness, section (4) studies the local and global stability, section (4) describes the MsDTM for finding the semi-analytical solutions and these results are discussed in section (6), and the last section contains a conclusion of this work.

Mathematical model
In this section, we introduce the mathematical model that describes the interactions between two types of predators (i.e Y and Z) and one type of prey (X). We aim to study the behavior of the animal population during the mating period under the following assumptions: • a male of the predators stays with the same female and they feed on the same prey during their mating period.
• the rate of growth for predators and prey are larger than the rate of decay by their natural death.
The modified diffusive predator-prey model in reference [22] studied the interaction of two predators (male and female from the same species) with one prey in the mating period. Herein, we extend this model to study two species of predators interacting with one prey during their mating period as follows where X(t), Y(t) and Z(t) denote densities of prey, predator type (I) and predator type (II), respectively. The terms XY and XZ refer to the interaction between prey and predators of type (I) and type (II) respectively. The terms XY 2 and XZ 2 represent the interaction between the male and female of predators and prey. All parameters are nonnegative and have the following interpretations • β 1 , β 8 , β 11 are the growth rates of X, Y, Z respectively.
• β 2 , β 9 , β 12 are the decay rates of X, Y, Z respectively due to natural death.
• β 3 , β 10 , β 13 are the decay rates of X, Y, Z respectively due to competition of food on the same species.
• β 4 , β 7 are the decay rates of prey due to the interaction between one predator and one prey.
• β 5 , β 6 are the decay rates of prey due to the interaction between two predators and one prey.

Positivity and boundedness
Theorem 1 All solutions of system (1) are nonnegative and bounded in some region O subject to nonnegative initial conditions in O.

Proof
Let assume the initial conditions at t = 0 are nonnegative and rewrite the system (1) as by integrating over [0, t], the following equalities are valid: and Hence, the solution (X(t),Y(t), Z(t)) is nonnegative for nonnegative initial data. Next, we prove that all solutions will remain bounded. We start by rewriting the system (1) as follows , (β 11 − β 12 )}, and γ 2 = min{β 3 , β 10 , β 13 }, thus, Hence, all the solutions of model (1) are nonnegative and bounded in the following region

Stability analysis
An equilibrium point of a dynamical system denotes a stationary condition for the dynamics. The local stability at equilibrium points indicates that the solutions remain near the equilibrium points with small perturbation, but can lead to a large change with large perturbation. On the other hand, the global stability means the solution remain nearby the equilibrium points with small or large perturbations. Thus, in this section, we will study the local and global stability of the equilibrium points.

Equilibrium points
To study the stability of the system (1), the equilibrium points are obtained by setting The solutions for the system (2) give the equilibrium points as follows: Note that the coexisting equilibrium point were excluded because it is negative points. Also, the trivial point was excluded because is not interesting in our investigation.  ; b 10 6 ¼ 0; b 13 6 ¼ 0;

Local stability
• (2): for the equilibrium point: ; b 10 6 ¼ 0; • (4): for the equilibrium point: The Jacobian matrix for the system (1) To study the stability of each equilibrium point, we substitute the point P i , i = 1, 2, 3, 4 into the matrix J and the result is J P . Then, looking for negative eigenvalues of J P implies the stability conditions for each equilibrium point.
(1): the equilibrium point: ð0; b 8 À b 9 b 10 ; b 11 À b 12 b 13 Þ is non negative if β 8 � β 9 , β 11 � β 12 , and the eigenvalue of J P 1 is negative if Þ is non-negative if β 8 � β 9 , β 11 � β 12 and the eigenvalues of J P 2 are nega- 0Þ is non-negative if β 8 � β 9 , β 11 � β 12 , and the eigenvalues of J P 3 are negative if b 1 � b 5 ðb 8 À b 9 Þ 2 þb 10 ðb 4 ðb 8 À b 9 Þþb 2 b 10 Þ b 2 10 ; b 10 6 ¼ 0 0Þ is non-negative if β 1 � β 2 , and the eigenvalues of J P 4 are negative if Thus, the points, ð0; b 8 À b 9 b 10 ; b 11 À b 12 b 13 Þ; ð0; 0; b 11 À b 12 b 13 Þ; ð0; b 8 À b 9 b 10 ; 0Þ and ð0; 0; b 11 À b 12 b 13 Þ are asymptotically stable associated to the presented stability conditions. The relevant ecological interpretation for stability analysis that the system has stable equilibrium point (1) if the growth rate of predators is greater than the natural death rate for both types of predators. The equilibrium point (2) is stable if the growth rate of the predator type (I) is less than the death rate and the growth rate of predator type (II) is greater than its death rate. In case of the opposite of this scenario, the equilibrium point (3) is stable. The equilibrium point (4) is stable if the growth rate of the prey is greater than the death rate which is a result of decreasing growth rate of predators type (I) and (II).

The multistage differential transform method
The idea of the multistage differential transform method is based on dividing the domain of time into subdomains [t i t i+1 ]. In each subdomain, we apply the differential transform method (DTM) where the initial condition at t i . The DT is defined as follows Let assume an ODE with initial condition at t 0 = 0 and a solution u(t). The operators in the ODE are transferred to DT operators as in Table 1. Then, we obtain the scheme in terms of U (k) and U(k + 1). The U(0) is the initial condition and U(k + 1) is obtained by applying U(k).
where N is the arbitrary iteration number. The readers are referred to the references [26] for more details.

Semi-analytical simulation
In this section, the proposed model will be computed by MsDTM. First step, we divide the , Z 0 (t i ) are the initial values and the X k+1 , Y k+1 , Z k+1 are obtained by applying DTM method in each subinterval as follows for k = 0, 1, . . ., α − 1. The solution of each case in Table 2 is plotted in Fig 1. The solutions clarify that if there are two types of predator such that male and female feed together on same prey during mating period, the predator who has low growth rate decays fast. However, the optimal ecological scenario is case (4) because it gives a periodic solution. The case (4) denotes that case when the growth rate of prey is larger than the growth rate of predators (I) and (II).

Conclusion
The article is a first contribution in the literature that produces two-predator-one-prey system with the effect of XY 2 term. It also solved the aforementioned mathematical model by MsDTM. The used method is fast, reliable and global convergent. The stability of the dynamic of the system was investigated. We found that the equilibrium points are asymptotically stable. The advantages of the used method have been proved in our previous work [26]. The method is more accurate than ADM, DTM, LADM and is faster than 4 th RK. In fact that the 4 th RK method discretizes the domain and its speed based on the step size of the discretization, but the speed of MsADM based on two factors: the length of subdomain and number of taken Taylor expansion terms. Therefore we can control the speed and accuracy of MsDTM by two factors.
The results proved that when two types of predators live in a same spot with the assumptions that male and female feed on the same prey during mating period, the ecological balance is achieved when the growth rate of prey is grater than both types of predators. Otherwise, the population of prey will decay fast and then the predator populations will decay as results of lack of food.
In future work, MsDTM can be used to solve different ODEs in the science and the presented model can be studied with fractional derivative.